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Abstract 

A detailed review of recent developments in Reynolds stress modeling for incompressible 
turbulent shear flows is provided. The mathematical foundations of both two-equation mod- 
els and full second-order closures are explored in depth. It is shown how these models can 
be systematically derived for two-dimensional mean turbulent flows that are close to equi- 
librium. A variety of examples are provided to demonstrate how well properly calibrated 
versions of these models perform for such flows. However, substantial problems remain 
for the description of more complex turbulent flows where there are large departures from 
equilibrium. Recent efforts to extend Reynolds stress models to non-equilibrium turbulent 
flows are discussed briefly along with the major modeling issues relevant to practical Naval 
Hydrodynamics applications. 
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1. Introduction 

Turbulent shear flows are of central importance for a variety of Naval Hydrodynamics 
applications ranging from flow around submerged bodies to free surface flows. Most of 
these turbulent flows are at extremely high Reynolds numbers - and in complex geometrical 
flow configurations - where the application of direct or large-eddy simulations are all but 
impossible for the foreseeable future. Reynolds stress models are likely to remain the only 
technologically feasible approach for the solution of these problems for the next few decades 
to come, if not beyond (see Speziale [1]). 

It is widely believed that Reynolds stress models are completely ad hoc, having no formal 
connection with solutions of the full Navier-Stokes equations for turbulent flows. While this 
belief is largely warranted for the older eddy viscosity models of turbulence, it constitutes a 
far too pessimistic assessment of the current generation of Reynolds stress closures. It will 
be shown how second-order closure models and two equation models with an anisotropic 
eddy viscosity can be systematically derived from the Navier-Stokes equations when one 
overriding assumption is made: the turbulence is locally homogeneous and in equilibrium. 
Moderate departures from equilibrium - where there are weak inhomogeneous effects - can 
then be accounted for in a relatively straightforward fashion. 

A brief review of zero and one equation models based on the Boussinesq eddy viscosity 
hypothesis will first be given in order to provide a perspective on the earlier approaches 
to Reynolds stress modeling. However, it will then be argued that since turbulent flows 
contain length and time scales which can change dramatically from one flow configuration 
to the next, two-equation models constitute the minimum level of closure that is physically 
acceptable. Typically, modeled transport equations are solved for the turbulent kinetic 
energy and dissipation rate from which the turbulent length and time scales are built up; 
this obviates the need to specify these scales in an ad hoc fashion for different flows. While 
two-equation models represent the minimum acceptable level of closure, second-order closure 
models constitute the highest level of closure that is currently feasible from a practical 
computational standpoint. It will be shown how the former models follow from the latter 
in the equilibrium limi t of homogeneous turbulence (see Speziale, Sarkar and Gatski [2] and 
Gatski and Speziale [3]). However, it will be demonstrated that the two-equation models 
which are formally consistent with second-order closures have an anisotropic eddy viscosity 
with strain- dependent coefficients - features that the most commonly used models do not 
possess. 

For turbulent flows that are only weakly inhomogeneous, full Reynolds stress closures can 
then be constructed by the addition of turbulent diffusion terms that are formally derived via 
a gradient transport hypothesis. Properly calibrated versions of these models are found to 
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yield a surprisingly good description of a wide range of two-dimensional mean turbulent flows 
that are near equilibrium. In particular, plane turbulent shear flows are accurately described 
with the stabilizing or destabilizing effect of a system rotation predicted in a manner that is 
quantitatively consistent with hydrodynamic stability theory. However, existing second-order 
closures are not currently capable of properly describing turbulent flows that are far from 
equilibrium and have major problems with wall-bounded turbulent flows. In regard to the 
latter point, it will be argued that we do not currently know how to properly integrate second- 
order closure models to a solid boundary with the no-slip condition applied. A variety of ad 
hoc wall damping functions are currently used that depend on the unit normal to (and/or 
the distance from) the wall - a feature that makes it virtually impossible to reliably apply 
these models in complex geometries. Consequently, in many applications of second-order 
closures to wall-bounded turbulence, the integration is carried out by matching to law of 
the wall boundary conditions, which do not formally apply to complex turbulent flows. The 
really disturbing feature here is that many of the commonly used second-order closures are 
not even capable of reproducing law of the wall results for an equilibrium turbulent boundary 
layer unless an ad hoc wall reflection term is added. This term typically depends inversely 
on the distance from the wall, further compromising the ability to apply these models in 
complex geometries. Entirely new approaches to the modeling of complex non-equilibrium 
and wall-bounded turbulent flows will be discussed briefly. 

A variety of illustrative examples involving turbulent shear flows will be provided in order 
to amplify the central points discussed in this paper. In addition, a special effort will be 
made to address the crucial issues in turbulence modeling that are relevant to practical Naval 
Hydrodynamics applications. 


2. Basic Equations of Turbulence 

We will consider the incompressible turbulent flow of a viscous fluid under isothermal 
conditions. The velocity field u,- and kinematic pressure P are solutions of the Navier-Stokes 
and continuity equations given by 

dP 


dvi dvi 
dt + Vj dx 


j dxi 
dvi 


+ vV Vi 


dxi 


= 0 


( 1 ) 


( 2 ) 


where v is the kinematic viscosity and the Einstein summation convention applies to repeated 
indices. As in all traditional studies of turbulence modeling, the velocity and kinematic 
pressure are decomposed into mean and fluctuating parts as follows: 

Vi = Vi + Ui, P — P + p (3) 
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where an overbar represents a Reynolds average. This Reynolds average can take a variety 
of forms for any flow variable </>: 


Homogeneous Turbulence 


Statistically Steady Turbulence 


6 — lim f d>(xA)d 3 x 

V— »oo VJv 

(f) = (j){t ) (Spatial Average) 

1 r T 

4 = <K x > t ) dt 

T—yoo 1 Jo 

(j) = <£(x) (Time Average) 


(4) 


(5) 


General Turbulence 


< l >= N m It 

AT— >oo Jy 


a=l 


( 6 ) 


(j) = <yfr(x,£) (Ensemble Average). 

In (6), ol represents a given realization of the turbulence. 

The Ergodic Hypothesis is assumed to apply. In a homogeneous turbulence, 

^ensemble — ^spatial 

whereas in a statistically steady turbulence, 

^ensemble — ^t ime* 


(7) 


( 8 ) 


For general turbulent flows that are neither statistically steady nor homogeneous, ensemble 
averages should be used (cf. Hinze [4] for a detailed discussion of these issues). 

The Reynolds- averaged Navier- Stokes and continuity equations take the form (cf. Hinze 

W) 


dvi _ dvi 


dt 


+ ^ = 


8P ,, drtj 

+ vV 2 Vi 


dxi 


dxj 


dvi 

dx: 


= 0 


where 


Tij — UiUj 


(9) 

( 10 ) 

( 11 ) 


is the Reynolds stress tensor. 
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The Reynolds-averaged Navier-Stokes equation is not closed until a model is provided 
that ties the Reynolds stress tensor to the global history of the mean velocity V{ in a 
physically consistent fashion. In mathematical terms, r,j is a functional of the global history 
of the mean velocity field, i.e., 

n-j(x,t) = ^i[v(x', *');*»*] 

x'eV, ' 

where Fij[ ‘ ] denotes a functional over space and time, and V represents the fluid volume. 
In (12), it is understood that there is an implicit dependence on the initial and boundary 
conditions for Ui and, hence, on those for the entire hierarchy of moments constructed from 
the fluctuating velocity. For the construction of Reynolds stress closures, it is typically 
assumed that the initial and boundary conditions for any turbulence correlations beyond the 
Reynolds stress tensor and dissipation rate merely serve to set the level of the length and 
time scales (see Lumley [5] and Speziale [1]). 


3. Zero and One Equation Models 


The Reynolds stress tensor can be decomposed into isotropic and deviatoric parts as 
follows: 

r.j = \k8„ + DTij (13) 

where the deviatoric part £>r,j is a symmetric and traceless tensor. Virtually all of the 
commonly used Reynolds stress models in this class are based on the Boussinesq hypothesis 
where it is assumed that 

D Tij = -v T 

given that vt is the eddy viscosity. For most incompressible turbulent flows, the isotropic 
part of the Reynolds stress tensor (§^) is not needed for the determination of the mean 
velocity field since it can simply be absorbed into the mean pressure P in (9). 

The eddy viscosity can be written as 


dvi dv; 
dxj dxi 


(14) 


= f (15) 

to 

where Iq is the turbulent length scale and to is the turbulent time scale — quantities that 
can vary dramatically with space and time for a given turbulent flow. In zero equation 
models , both £ 0 and t 0 are specified algebraically by empirical means. The first successful 
zero equation model based on the Boussinesq eddy viscosity hypothesis was Prandtl’s mixing 
length theory (see Prandtl [6]). In Prandtl’s mixing length theory, 


VT = ll\ 



(16) 
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where £q — ny is the mixing length, k is the von Karman constant, and y is the normal 
distance from a solid boundary. This representation is only formally valid for thin turbulent 
shear flows - near a wall - where the mean velocity is of the simple unidirectional form 
v = u(y) i. 


Several decades later, this simple mixing length model was generalized to multi-dimensional 
turbulent flows. Three alternative tensorially invariant forms have been proposed: 

Smagorinsky [7] Model 

v T = ll{2 Susa) 1 !* 

(17) 

Cebeci- Smith [8] Model 

2 / dvi dvi \ 1/2 
T 0 \dx j dx j J 

(18) 

Baldwin- Lomax [9] Model 

v T = £ 2 q (a\-o7i) 1/2 

(19) 


where Sij = ~{dvi/dxj + dvj/dxi ) is the mean rate of strain tensor and a? = V x v is the 
mean vorticity vector. The former model has been primarily used as a subgrid scale model for 
large-eddy simulations whereas the latter two models have been used for Reynolds- averaged 
Navier- Stokes computations in aerodynamics. Each of these models reduces to the simple 
mixing length formula (16) in the thin shear flow limit. However, they suffer from the same 
deficiency as the original mixing length model in their need for an ad hoc specification of 
the turbulent length scale £ 0 - a task that is all but impossible to do reliably in complex 
turbulent flows. 

Beyond the length scale specification problem with zero equation models, there is another 
criticism that can be raised: it is not physically consistent to build up the turbulent velocity 
scale from the mean velocity gradients as done in (17) - (19)- The proper measure of the 
turbulent velocity scale is the intensity of the turbulent fluctuations (i.e., we should take 
v 0 = K 1 / 2 where v 0 = £o/to). Hence, a more physically consistent representation for the 
eddy viscosity is given by 

ur = K'/Ho. ( 20 ) 

Prandtl [10] - who expanded on many of the earlier ideas of Kolmogorov [11] - developed a 
one-equation model based on (20) wherein a modeled transport equation for the turbulent 
kinetic energy was solved. Subsequent to this early work, a variety of researchers have 
proposed one-equation models along these lines for near-wall turbulent flows (cf. Norris and 
Reynolds [12] and Rodi, Mansour and Michelassi [13]). 
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One-equation models based on the solution of a modeled transport equation for the 
turbulent kinetic energy still suffer from one of the major deficiencies of mixing length models: 
they require the ad hoc specification of the turbulent length scale which is virtually impossible 
to do reliably in complex three-dimensional turbulent flows. Recently, one-equation models 
have been proposed based on the solution of a modeled transport equation for the eddy 
viscosity v t (see Baldwin and Barth [14] and Spalart and Allmaras [15]). These models do 
alleviate the problem of having to specify the turbulent length scale in the definition of the 
eddy viscosity (20). Nonetheless, an ad hoc specification of length scale must be made in the 
destruction term within the modeled transport equation for v? which depends empirically 
on the distance from the wall. 

This leads us to one of the central points of this paper: the turbulent length and time 
scales (4, t 0 ) are not universal; they depend strongly on the flow configuration under consid- 
eration. Consequently, two-equation models - wherein transport equations are solved for two 
independent quantities that are directly related to the turbulent length and time scales - 
represent the minimum acceptable level of closure. In the most common approach, the turbu- 
lent length and time scales are built up from the turbulent kinetic energy K and dissipation 
rate e (i.e., £ 0 oc K 3 ^ 2 /e, t 0 oc K/e) with modeled transport equations solved for K and e. 
These two-equation models should be formulated with a properly invariant anisotropic eddy 
viscosity that is nonlinear in the mean velocity gradients. The standard Boussinesq eddy 
viscosity hypothesis makes it impossible to properly describe turbulent flows with: (a) body 
force effects arising from a system rotation or from streamline curvature, and (b) flow struc- 
tures generated by normal Reynolds stress anisotropies (e.g., secondary flows in non-circular 
ducts). 

At this point, it would be useful to comment on the most sophisticated level of Reynolds 
stress closure that is now practical. Limitations in computer capacity, and issues of numerical 
stiffness, appear to make second-order closure models — wherein modeled transport equations 
are solved for the individual components of the Reynolds stress tensor along with a scale 
equation - the highest level of closure that is currently feasible for practical computations. 


4. Turbulent Transport Equations 

The transport equation for the fluctuating velocity u,-, which is obtained by subtracting 
(9) from (1), takes the form 


dui _ dui 
+ Vj 


dt 


dx , 


= —Ui 


dui 


3 dxj 
+ vV 2 Ui - 4 * 


u 


dvi 


3 

dra 


dxj 


dp 

dx{ 




( 21 ) 


dxj' 
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This equation can be written in operator form as 

Mui = 0. 


( 22 ) 


The Reynolds stress transport equation is obtained by constructing the second moment 


UiMuj + UjAfui = 0 . 


(23) 


Its full form is given by (cf. Hinze [4]) 


drij _ dTij dvj &Vi 

at +Vk dx k ~ Tik dx k Tik dx k + * iJ 


dCijk __2 
-e.'j - -5— + "V Ti y 


(24) 


dx k 


In (24), 


$ij = p 


du; duj 

+ 3 


dxj dx{ 

_ dui duj 

£ ' 3 = 2v d^ k dx\ 


(25) 

(26) 

Cijk — 'U/i'lLjlLk “f* JpXLiSjfc pUj&ik (^) 

are, respectively, the pressure-strain correlation, the dissipation rate tensor and the turbulent 
diffusion correlation. 

The transport equation for the turbulent kinetic energy K = is obtained by con- 
tracting (24): 

+ '4f- 7 ’ — s-g«**;+w) + - v ^ <*> 


dK 


dt 


where 


V = ~T t 


dv i 


e = v 


tJ dxj 

dui dui 
dxj dxj 


(29) 


(30) 


are, respectively, the turbulence production and the turbulent dissipation rate. By con- 
structing the moment 

(31) 


2 = 0 , 


dxj dxj 

a transport equation for the turbulent dissipation rate e can be obtained. This equation 
takes the form [1] 

% + v~ = V c -$ c +V s + vV 2 e (32) 

dt oxi 
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where 


—2v 


-2u 


duk duk dvi n dui duj dvi 


dxi dxj dxj dxk dx k dxj 

(33) 

duk du k dui o dui d' l Vi 

dxi dxj dxj k dxj dxjdxk 


2 d 2 ui d 2 Ui 

— Zu 

dx jdiisk dtBjdtCfz 

(34) 

( dp duj \ d 1 dui dui \ 

\dxidxi) V dxj \ 3 dx k dxk) 

(35) 


are, respectively, the production, destruction and turbulent diffusion of dissipation. 

Both two-equation models and second-order closure models are obtained by modeling 
the Reynolds stress transport equation (24) and the dissipation rate transport equation 
(32). Second-order closures are obtained by modeling the full Reynolds stress transport 
equation. Two-equation models are formally obtained by assuming that the turbulence is 
locally homogeneous and in equilibrium; the Reynolds stress anisotropies are then derived 
algebraically from (24) and a modeled version of (28) for the turbulent kinetic energy is 
solved. 


5. Two-Equation Models 

It will now be shown how two-equation models can be systematically derived from the 
Reynolds stress transport equation. As alluded to earlier, two-equation models — with an 
algebraic representation for the Reynolds stresses — are obtained by assuming that the tur- 
bulence is locally homogeneous and in equilibrium. Hence, we start with the Reynolds stress 
transport equation for homogeneous turbulence given by: 


T ij — 


-T ik 


dvj 

dx k 


Ovi 

Tik dx k + 


6ij . 


Since, the fluctuating pressure p is a solution of the Poisson equation 

V 2 p = 


dui duj dvi duj 


dxj dxi dxj dxi 
it follows that the pressure-strain correlation can be written in the form 

dvk 


4 ? ij — Mijki 


dxi 


In (38). 


*-£/// id- 


du% duf ( dui ( duj 


+ 


x*| dxf dxi \dxj dxi 


d 3 


x * 


(36) 


(37) 


(38) 


(39) 
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( 40 ) 


are, respectively, the slow and rapid terms which are obtained by implementing the Green’s 
function solution of (37) for an infinite flow domain. 

In the developments to follow, extensive use will be made of the Reynolds stress and 
dissipation rate anisotropy tensors defined as 


_r,,-§ KSa 

°tj — 


2 K 


e 1 c 

c, tJ 

d ' j = 2e - 


(41) 


(42) 


respectively (see Lumley [16] and Reynolds [17]). Furthermore, use will be made of the trans- 
port equation for the turbulent kinetic energy which is exact for homogeneous turbulence: 


K = P — e 


(43) 


Eq. (43) is obtained by contracting (36). The direct substitution of (38) - (42) into (36) 
yields the Reynolds stress transport equation 


b„ = -~s„ - (b ik ^ + b, k ^ - ^b k r~s 


2U 
3 


dx k 


dvi 

Jjk 8^ k 


2. dv t 

3'~ 9*r ,j 


(44) 


^ 2 K-M-ijki 


dvk 

dxi 


—dij 

K J 


given in terms of the anisotropy tensors alone. In (44), Aij = A ij/e and Mijki = Mijki/K 
are the dimensionless slow and rapid pressure-strain terms. 

The fundamental assumptions underlying two-equation models are that the turbulence 
is locally homogeneous and an equilibrium state is reached where 

K 

bij ) d^ , Aij , Mijki j ~ 

attain constant values that are largely independent of the initial conditions. In general, Aij 
and Mijki are functionals, in wavevector space k, of the energy spectrum tensor i£ tJ f(k,t) 
where ^ 

r,, = JJ J E„(k,t)d 3 k (45) 

— oo 

(cf. Reynolds [17]). This has prompted turbulence modelers to construct one-point models 
for Aij and Mijki of the form (Lumley [16]) 


Aij — Aij ( b ) , Ad ijki Ad ijki ( b ) . 


(46) 
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It should be said at the outset that models of the form (46) cannot be expected to apply to 
general homogeneous turbulent flows since nonlocal effects in wavevector space are neglected; 
it is well known that Mijki is of the form (cf. Reynolds [17]) 


M i]kl ~ JJ I ^E kl (k,t)d 3 k. 


(47) 


However, for homogeneous turbulent flows that are in equilibrium, there is evidence to 
suggest that Aij,M,ijki and b{j achieve constant values that are independent of the initial 
conditions as alluded to earlier (homogeneous shear flow represents a prime example; see 
Tavoularis and Corrsin [18]). Any constant tensor can be written as a finite expansion in 
three linearly independent vectors that are also constant. Since bij is a symmetric tensor, 
its eigenvectors are linearly independent; hence, (46) is expected to be formally valid for a 
homogeneous turbulence that achieves this type of structural equilibrium. 

Speziale, Sarkar and Gatski [2] showed that for two-dimensional mean turbulent flows 
that are homogeneous and in equilibrium, the pressure-strain correlation reduces to the 
simple general form: 

dv k 


*ij = eAij{b) 4 - KM ijU (b) 


dx { 


CiSbij ~f* C/-y£ ^bikbkj 

-\-G2KSij -f C3K (bikSjk + bjkSik 

2 _ \ 

-bkiSkiSij) C 4 K (bikSjk 4 * bjktOik) 


(48) 


where 


S *=2 


U) 


IJ 


' dvi dvj 
dxj + 

dvi 
dxj 



(49) 


(50) 


are, respectively, the mean rate of strain tensor and the mean vorticity tensor. In (48), 
C\ — C 4 are constants that are not necessarily universal; in principal, their specific numerical 
values can vary from one flow to the next. However, it is encouraging to note that, consis- 
tent with its definition (47), the basis expansion (48) has a rapid part that is linear in T ij 
and, hence, linear in the energy spectrum tensor. It is only in the limit of two-dimensional 
mean turbulent flows that the general basis expansion for ( 46)2 satisfies this linear con- 
sistency condition - a result of the fact that II, III and fc 33 achieve universal equilibrium 
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values in the two-dimensional limit (Speziale, Sarkar and Gatski [2]). For uniformly strained 
turbulent flows near equilibrium, there is substantial evidence from physical and numerical 
experiments to suggest that the quadratic return term in (48) (with coefficient <7*) can be 
neglected without introducing an appreciable error. Then, the representation (48) becomes 
completely linear in the Reynolds stress tensor. This allows for the superposition of solutions 
and m aint ain s consistency with the linearity of the rapid pressure-strain correlation in the 
energy spectrum tensor — a property that follows from its definition as stated above. In the 
opinion of the author, this constitutes the primary reason for the relative success that (46) 
has had in the description of two-dimensional mean turbulent flows that are near equilib- 
rium. The applicability of (46) to non-equilibrium turbulent flows or to three-dimensional 
mean turbulent flows is highly debatable. In regard to the latter case, the general basis rep- 
resentation for (46)2 is highly nonlinear in 6 tJ - (see Lumley [16], Reynolds [17] and Speziale 
[1]) — and, therefore, nonlinear in the energy spectrum tensor - in violation of (47). 

If we neglect the anisotropy of dissipation, then in the equilibrium limit where 6 tJ - = 0, Eq. 
(44) reduces to the following linear system of algebraic equations (see Gatski and Speziale 

[ 3 ])= 


b*j — -S i:i - b* ik S jk - b* h S ik + ~b* kl S kl 8ij 
+b* ik W kj + b* jk W ki 


where 


S*,- = \gj { 2 - C,)S tl 
w *ij = 1 * 7(2 - Ctpa 


*«= 


C3-2 

C, -i 


^ij 


(Oi.r A " 1 

9 -{y + 7~ 1 ) ' 

For turbulent flows in non-inertial frames of reference, Coriolis terms must be added to the 
right-hand-side of (44) along with a non-inertial correction to the pressure-strain correlation 
model (48). As shown by Gatski and Speziale [3], these terms can be accounted for exactly 
by simply replacing (53) with the extended expression 


(51) 

(52) 

(53) 

(54) 

(55) 


Wi = \g*(2 - ft) 


(jJij -f 


(C A - 

\c 4 - 


, ..o 


(56) 


where e m ji is the permutation tensor and O m is the angular velocity of the reference frame 
(in an inertial frame of reference, where f2 m = 0, the expression (56) reduces to (53)). 
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Equation (51) constitutes a set of linear algebraic equations for the determination of b *- 
in terms of ~S*j and W*j; the solution to (51) is of the general mathematical form 


b* =f{T,W). 


(57) 


As first suggested by Pope [19], the general solution to the implicit algebraic stress equation 
(51) is of the form: 

10 

b* = £ G W T (a) (58) 


where 


A=1 


= gp, 


T< 6 > = W* 2 S*+S*W 


*■=■==* 2 - 


— ?{s w }i 


T< 2 > = S*W*- W*S*, T< 7 > = W*S*W* 2 

-W* 2 S*W* 


T( 3 > = s ' 2 - 1{S* 2 }I, T< 8 > = s*W*s * 2 


-s w s 


(59) 


X (4) = w* 2 - §{W* 2 }I, T( 9 ) = W* 2 S* 2 +S* 2 W 2 

-§{s’* 2 w* 2 }i 

T (5) _ w*S* 2 - S* 2 W*, T( 10 > = W*s* 2 W * 2 

-W* 2 S* 2 W* 

are the integrity bases ({•}■ denotes the trace). Pope [19] only obtained the solution to (51) 
corresponding to the Launder, Reece and Rodi [20] model simplified to two-dimensional mean 
turbulent Hows in an inertial frame — a case for which the calculations become much simpler 
since only the integrity bases are linearly independent. Gatski and Speziale [3] 

showed that the general solution (58) for three-dimensional turbulent flows is as follows: 


G« = -1(6 - 3^ - 21t7 2 
—2773 + 30774)//), 

G< 2 > = -(3 + 3771 - 6772 
+2773 + 6774)/!), 

G< 3 > = (6 - 377! - 127)2 
-2773 - 6774)/!), 

G< 4 > = —3(377! + 2773 

+6774)/!), 

G< 5 > = -9/D, 


G< 6 > = -9/D 
G( 7 > = 9/D 
G< 8 > = 9/D 
G< 9 > = 18/D 
G< 10 > = 0 


( 60 ) 


7 , 15 „ 2 , 2 

D = 3 - -7)! + 7)! 2 - y7)2 - Sl)lV2 + 3t) 2 - 7) 3 + -T)^ 
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(61) 

(62) 


-2t) 2 7}3 4 217/4 4 247/5 4 27/47/4 - 67/ 2 7/ 4 

t/ 4 = {r 2 >, t/ 2 = {w* 2 }, 7/3 = {r 3 >, 

7/4 = {S*r 2 }, 7/s = {S* 2 W* 2 }. 

While the results provided in (60) - (62) constitute the general solution of (51) for three- 
dimensional turbulent flows, questions can be raised about its overall usefulness. As alluded 
to earlier, (51) is based on the use of (48) which is only formally valid for two-dimensional 
mean turbulent flows that are near equilibrium. For two-dimensional mean turbulent flows, 
(60) - (62) simplifies substantially to the form 


6?. 

*7 


3 - 2t/ 2 4 6f 2 


& + r ik w kj + 


-2 ( s* k slj - 


(63) 


where 

v = (KK) 1/2 > ( = (Knr /2 (64) 


By making use of (52) - (55), we can write (63) in terms of the Reynolds stress tensor as 
follows: 


T,J 3 KS,j 3-2ij 2 + 6£ 2 


a i 






T(3 __ __ 

+a2 nr (S'i*W%- + SjkWti) 


(65) 


where 



(SikSkj — -SkiSuSij) 


= g (| - C 2 ) (66) 

“2 = ~g 2 — C 2 ) ( 2 — C 4 ) (67) 

<*3 = g 2 (| - C,) (2 - ft) (68) 

The coefficients a 1? a 2 and a 3 are not constants but rather are related to the coefficients 
C\ — C 4 and g. In mathematical terms, they are “projections” of the fixed points of Aij and 
Mijki onto the fixed points of , which can vary from one flow to the next. However, for 
2 -D turbulent flows, C\ — C 4 can be approximated by constants due to the linear dependence 
on b{j which allows us to use superposition . 

Gatski and Speziale [3] evaluated C\ — C 4 using the SSG second-order closure which will 
be discussed later; this model was calibrated largely based on the use of data for homogeneous 
shear flow (see Table 1). 
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Equilibrium 

Values 

LRR 

Model 

SSG 

Model 

Data 

(^11)00 

0.158 

0.204 

0.201 

(^22)00 

-0.123 

-0.148 

-0.147 

(^12)00 

-0.187 

-0.156 

-0.150 

(SK/e) x 

5.32 

5.98 

6.08 


Table 1. Comparison of the predictions of the Launder, Reece and Rodi (LRR) model and 
the Speziale, Sarkar and Gatski (SSG) model with the experimental data of Tavoularis and 
Corrsin [18] for homogeneous turbulent shear flow. 

The constant values that are taken for C\ — C 4 are given by (see Gatski and Speziale [3]): 

C x = 6.80, C 2 = 0.36, G 3 = 1.25, C 4 = 0.40 (69) 

It should be noted at this point that if (63) is applied to turbulent flows that are far from 
equilibrium, singularities can arise through the vanishing of the denominator containing 7 7 
and £ (it is straightforward to show that this cannot happen in equilibrium turbulent flows). 
Hence, this model needs to be regularized before it is applied to complex turbulent flows that 
are not in equilibrium. This can be accomplished via a Pade type approximation whereby 

3 „ 3(1 +t? 2 ) 

3 - 2t7 2 + 6£ 2 ~ 3 + 7 ] 2 + 6f V + 6£ 2 

(see Gatski and Speziale [3]). It is a simple matter to show that (70) constitutes an ex- 
cellent approximation for turbulent flows that are near equilibrium and, unlike the original 
expression, is a bounded and non-negative function for all values of 7 ] and £. 

The representation (63) constitutes an anisotropic eddy viscosity model of the general 
form 

b. } = a tjkl g (71) 

where the fourth rank tensor a ijki depends on the symmetric and antisymmetric parts of 
the mean velocity gradients. Quadratic models of this type have recently been obtained 
by Yoshizawa [21], Speziale [22] and Rubinstein and Barton [23] based on two-scale DIA, 
continuum mechanics and RNG based techniques, respectively (in regard to the latter, see 
Yakhot and Orszag [24]). Furthermore, it must be noted that while the traditional implicit 
algebraic stress models such as that due to Rodi [25] (which is of the general form (51)) 
have an explicit solution of the form (63), they are ill-behaved and can give rise to diver- 
gent solutions when applied to non-equilibrium turbulent flows. This explains why previous 
anisotropic corrections to eddy viscosity models have only had limited success: 
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(i) A quadratic expansion is not adequate; the coefficients should depend nonlinearly on 
rotational and irrotational strain rates. 

(ii) Only the regularized explicit solution to algebraic stress models - which has just re- 
cently emerged - has the proper such dependence. Traditional algebraic stress models 
are ill-behaved and should not be applied to complex turbulent flows that are signifi- 
cantly out of equilibrium. 

If we have a clear cut separation of scales where 

* 7 , t < 1 


then (65) reduces to the eddy viscosity model 

m = \KSii - 2C„— Sij (72) 

which forms the basis for the standard K — e model of Launder and Spalding [26]. However, 
in basic turbulent shear flows, we do not have a separation of scales: 77 and £ are of order one. 
Nonetheless, there are some circumstances where (65) yields results that are comparable to 
(72). For example, in the logarithmic region of an equilibrium turbulent boundary layer, the 
explicit algebraic stress model (65) yields 


' xy 


M e dy 


for the shear stress, where 


C„ n 0.094 


(73) 

(74) 


given that dvi/dxj = du/dy Sn 6 j 2 - This is virtually identical to the standard K — e model 
which, for this case, yields (73) with C ^ — 0.09. Of course, for more complex turbulent 
flows the models are substantially different; unlike the standard K — e model, the explicit 
algebraic stress model has a strain-dependent eddy viscosity and anisotropic eddy viscosity 
terms. 

In order to achieve closure, a modeled transport equation for the turbulent dissipation 
rate e is needed. For homogeneous turbulence, the exact transport equation (32) for the 
turbulent dissipation rate reduces to: 

dvi ( c ) dv t o du k du k dui 

e ~ ~ £ii dxj £ij dxj dxidxjdxj 

d 2 Uj d 2 Uj 
dxjdx k dxjdx k 
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where £ tJ is the turbulent dissipation rate tensor defined in (26) and 


.(<0 _ 


'ij 


2v 


duk duk 


( 76 ) 


dxi dxj 

is the complementary dissipation rate tensor. If we introduce the anisotropy of dissipation 
tensors 

2 

(77) 


dij — 


& c 

ij 


(c) _ Cjj 


g *’i 

2e 

(c) - ic& 




2e 


(78) 


(where e = | ea = |e£^), a simple closure can be developed for the production of dissipation 
terms in (75). Here, it is assumed that 

dij — C d bij , 


4 ? = c. 


d bij j 


(79) 


which physically implies that the anisotropy of dissipation is proportional to the anisotropy 
of the Reynolds stresses due to the fact that the former follows from the latter as a result 
of the energy cascade from large to small scales. Results from Direct Numerical Simulations 
(DNS) of homogeneous shear flow (Rogers, Moin and Reynolds [27]) only provide justification 
for (79) as, at best, a low order approximation. 

The third correlation on the right-hand-side of (75) can be written in the form 

.2 


2v 


where 


duk duk dui __ 7 

dx i dxj dxj 3\/i5 

6i/l5 


c Rl/2- 


Sk — 


duk duk duj 
dsci ()xj dx j 


( du m 8u r , 

dx jx dxj 


3/2 


(80) 


(81) 


is the generalized velocity derivative skewness and Rt = K 2 /ve is the turbulence Reynolds 
number. For isotropic turbulence, (81) reduces to the classical definition of the velocity 
derivative skewness which is given by Sk — — (du/dx) 3 /[(du/dx) 2 } 3 / 2 (here we define the 
skewness with the negative gauge). In spectral space, the destruction of dissipation term on 
the right-hand-side of (75) behaves as follows: 


2v* 


d 2 Ui d 2 Ui 


poo 

2i/ 2 / k 4 E(k,t)dk 
Jo 


(82) 

dxjdxk dxjdxk Jo v 7 ' 

where E(k,t) is the three-dimensional energy spectrum. Consequently, most of the contri- 
butions to this term occur at high wavenumbers where the energy spectrum scales with the 
Kolmogorov length scale, Ik = i/ 3 / 4 /e 1//4 . With this Kolmogorov scaling, it follows that 

.2 


2v l 


d 2 uj d 2 uj _ 7 

dxjdxk dxjdxk Zy/YS 


GkR^J 2 ^ + G £ 2 


(83) 
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(see Speziale and Bernard [28]). The direct substitution of (79), (80) and (83) into (75) 
yields the transport equation 


i = ~C e i-rr 


dvi 


K t3 dxj 


+ s7i5 (5v ~~ Gk ^‘ /2 W - C *K 


(84) 


where C £ \ — Cd + C% (in general homogeneous turbulence, C £ i, (7 e 2, Sk, and Gk can be 
functions of time). For equilibrium turbulent flows at high Reynolds numbers, 


Sk = Gk 


(85) 


and C £ i and C £ 2 can be approximated as constants (when (85) is not valid, then e changes 
on the Kolmogorov time scale - an extremely rapid change at high Reynolds numbers that 
constitutes a non- equilibrium flow situation). This leads us to the commonly used modeled 
dissipation rate equation for homogeneous turbulence: 

■ _ _ r e 9vi p f! 

Ccl K ,J dx } c2 K 

with C £ 1 and C £ 2 taken to be constants. Typically, C £ 2 is determined from the decay of 
isotropic turbulence; for isotropic decay, (86) implies that (cf. Speziale [1]) 

K ~ . (87) 

The most cited experimental data [29] indicates that the exponent of the decay law (87) has 
a mean value of approximately 1.2; this implies a value for C £ 2 ~ 1.83. In practice, a value 
of C e2 — 1.90 has been more commonly used starting with Launder, Reece and Rodi [20]. 
This typically has been used with a value of C £ 1 = 1.44 based on a calibration with a range 
of benchmark turbulent shear flows. 

Recently, Speziale and Gatski [30] showed that when the effects of anisotropic dissipation 
are more rigorously accounted for, a variable C £ \ results that is of the form C £ \ = G e i(77,£). 
This form is obtained by starting with a modeled transport equation for the full tensor dissi- 
pation £ij. An algebraic equation - analogous to that obtained from the ASM approximation 
for the Reynolds stress - is arrived at when the standard equilibrium hypothesis 

dij = 0 
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is invoked. For two-dimensional mean turbulent flows, it has tbe exact form: 


dij = -2 C, 


lie 






11 


,J yCes + V/e- 1 

7 30 ^ 2 \ 

11 n -~) (W 


( 88 ) 


+ 


c £5 -\-V js 


~S S^Sij] 

^ mn mn J 


where 


Cf “ 15(C a + V/e - 1 ) 


[1 




Ti a 3 4 - ^ 


li 


\C e s + P/e-l. 


_ JL 

li 3 li 


3 \c e5 + 7>/ £ -i; 


-l 


^ JiT 






Sij — Sij ^ , uj- — a;,- 

and (7 e5 and a 3 are constants (Speziale and Gatski [30]). The substitution of these algebraic 
equations into the contraction of the £ tJ - transport equation yields the scalar dissipation rate 
equation (86) with 


Cei — 1 + 


i5a, 


(1 + a)(a 5 + Oi 7 / 2 — i) 


(89) 


where 


(Ce5 + Crf - l) 2 + Pie - \PW 

V = (25* i 5* i ) 1/2 , C = (2 w| 3 w*,') 1/2 

3 (U 16\ 7 i 

" = 4 VlT a3 _ 33/ ’ ^-H Q3 + IT 

15 1 

P 2 = — , as ~ 5, a 3 « 0.6. 

The constants a 3 and as were evaluated using DNS results for homogeneous shear flow 
(Rogers, Moin and Reynolds [27]). 

For two-dimensional turbulent shear flows that are in equilibrium, (89) yields 


a! ~ i-4 
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which is remarkably close to the traditionally chosen constant value of C £ \ — 1.44. It is 
interesting to note that an alternative variable C £ i of the form C £ \ = C £ i(tj) was recently 
proposed by Yakhot et al. [31] based on a heuristic Pade approximation. However, the model 
of Speziale and Gatski [30] depends on rotational as well as irrotational strain rates ( 77 ,^). It 
has long been recognized that the dissipation rate is dramatically altered by rotations. The 
results of Speziale and Gatski [30] clearly show that this effect can be rationally incorporated 
by accounting for anisotropic dissipation. To the best knowledge of the author, this model 
constitutes the first systematic introduction of rotational effects into the scalar dissipation 
rate equation. Previous attempts to account for this effect (see Raj [32]; Hanjalic and 
Launder [33]; and Bardina, Ferziger and Rogallo [34]) were largely ad hoc. 

For weakly inhomogeneous turbulent flows that are near equilibrium, we can extend the 
K and e transport equations by the addition of turbulent diffusion terms that are obtained 
by a formal expansion technique: 


dK __ T - d 

— - + v • VK = V - e + 
at oxi 

de _ _ e e 2 

— + v - Ve = C el —V - C,_ 2 — + 



+ vV 2 K 


+ vV 2 e 


(90) 


(91) 


<t £ dx i 

where < 7 * and a £ are constants that typically assume the values of 1.0 and 1.3, respectively. 

This model can be integrated directly to a solid boundary, where the no-slip condition is 
applied, without the need for ad hoc wall damping functions. It is only necessary to remove 
the singularity in the destruction of dissipation term 


a 


e2 


K 


on the right-hand-side of (91). Durbin [35] argued that this expression should be replaced 
with the term 



where T is the turbulent time scale. For high Reynolds number turbulence, T = Kfe\ for 
low Reynolds number turbulence near a wall, the turbulent time scale is proportional to the 
Kolmogorov time scale, i.e., T oc \JvJs- These considerations lead Durbin [35] to propose 
the expression 


T = max 


K 

. e 


)C K \jvle 


where Ck is a constant of order one. A damping function, however, can also be used. 
Namely, we can take the destruction term to be 


-C £2 f 2 


K 
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where / 2 is a wall damping function which, for example, can be chosen to be of the form 

/ 2 = 1 - exp(— i? y /10) 

where R y = K^^y/v is the turbulence Reynolds number based on the distance y from the 
wall. No wall damping is needed in the eddy viscosity; the strain- dependent terms in the 
eddy viscosity provide natural damping as the wall is approached (see Speziale and Abid 
[36]). 

We will now consider several non-trivial applications of the two-equation model discussed 
herein which can be referred to as an explicit algebraic stress model (ASM) based on the SSG 
second-order closure. The first case that will be considered is homogeneous shear flow in a 
rotating frame (see Figure 1). In this flow, an initially isotropic turbulence (with turbulent 
kinetic energy Kq and turbulent dissipation rate £o) is suddenly subjected to a uniform shear 
with constant shear rate >9 in a reference frame rotating steadily with angular velocity Cl. In 
Figures 2(a)-2(c), the time evolution of the turbulent kinetic energy predicted by this new 
two-equation model is compared with the large-eddy simulations (LES) of Bardina, Ferziger 
and Reynolds [37], as well as with the predictions of the standard K — e model and the full 
SSG second-order closure. From these results, it is clear that the new two-equation model 
yields the correct growth rate for pure shear flow ( Cl/S = 0) and properly responds to the 
stabilizing effect of the rotations Cl/S = 0.5 and Cl/S = —0.5. These results are remarkably 
close to those obtained from the full SSG second-order closure as shown in Figure 2. In 
contrast to these results, the standard K — e model overpredicts the growth rate of the 
turbulent kinetic energy in pure shear flow ( Cl/S = 0) and fails to predict the stabilizing 
effect of the rotations illustrated in Figures 2(b)-2(c). Since the standard K - e model 
makes use of the Boussinesq eddy viscosity hypothesis, it is oblivious to the application 
of a system rotation (i.e., it yields the same solution for all values of Cl/S). The new 
two-equation model predicts unstable flow only for the intermediate band of rotation rates 
— 0.09 < Cl/ S < 0.53; this is generally consistent with linear stability theory that predicts 
unstable flow for 0 < Cl/S < 0.5. 

In Figure 3, the prediction of this new two-equation model for the mean velocity profile 
in rotating channel flow is compared with the experimental data of Johnston, Halleen and 
Lezius [38] for a rotation number Ro — 0.068. It is clear from these results that the model 
correctly predicts that the mean velocity profile is asymmetric in line with the experimental 
data - an effect that arises from Coriolis forces. In contrast to these results, the standard 
K — e model incorrectly predicts a symmetric mean velocity profile identical to that obtained 
in an inertial frame (the standard K — e model is oblivious to rotations of the reference frame, 
as alluded to above). As demonstrated by Gatski and Speziale [3], the results obtained in 
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Figure 3 with this new two-equation model are virtually as good as those obtained from a full 
second-order closure. This is due to the fact that a representation is used for the Reynolds 
stress tensor that is formally derived from a second-order closure (the SSG model) in the 
equilibrium limit. It is now clear that previous claims that two-equation models cannot 
systematically account for rotational effects were erroneous. 

Two examples will now be presented that illustrate the enhanced predictions that are ob- 
tained for turbulent flows exhibiting effects arising from normal Reynolds stress differences. 
Here, we will show results obtained from the nonlinear K — e model of Speziale [22]. For 
turbulent shear flows that are predominantly unidirectional, with secondary flows or recircu- 
lation zones driven by small normal Reynolds stress differences, a quadratic approximation of 
the anisotropic eddy viscosity model discussed herein collapses to the nonlinear K — £ model 
(see Gatski and Speziale [3]). In, Figure 4, it is demonstrated that the nonlinear K — e model 
predicts an eight-vortex secondary flow, in a square duct, in line with experimental obser- 
vations; on the other hand, the standard K — £ model erroneously predicts that there is no 
secondary flow. In order to be able to predict secondary flows in non-circular ducts, the axial 
mean velocity v z must give rise to a non-zero normal Reynolds stress difference r yy — r xx (see 
Speziale and Ngo [39]). This requires an anisotropic eddy viscosity (any isotropic eddy vis- 
cosity, including that used in the standard K — £ model, yields a vanishing normal Reynolds 
stress difference which makes it impossible to describe these secondary flows). 

In Figure 5, results obtained from the nonlinear K — £ model are compared with the 
experimental data of Kim, Kline and Johnston [40] and Eaton and Johnston [41] for turbulent 
flow past a backward facing step. It is clear that these results are excellent: reattachment 
is predicted at x/H ~ 7.0 in close agreement with the experimental data. In contrast to 
these results, the standard K — £ model predicts reattachment a t x/H « 6.25 - an 11% 
underprediction. This error predominantly results from the inaccurate prediction of normal 
Reynolds stress anisotropies in the recirculation zone as discussed by Speziale and Ngo 
[39] . As alluded to above, the new two-equation model can be integrated directly to a solid 
boundary with no wall damping. In Figure 6, the skin friction coefficient obtained from this 
model - plotted as function of the Reynolds number based on the momentum thickness, Re 
— is compared with experimental data and with results obtained from the K — £ model with 
wall damping. Clearly, the results are extremely good. 
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6. Second-Order Closure Models 

These more complex closures are based on the full Reynolds stress transport equation 
with turbulent diffusion: 


dTij _ dTij 

dt k dx k 


dvj _ dvi x ^ 

Tik ^ Tjk „ T *&ij 


dxi 


dx k 


2 C dCijk yy2 

-DSij - -e8ij - -f uV Tij 


(92) 


where D^ij is the deviatoric part of the dissipation rate tensor. Full second-order closure 
models are needed for turbulent flows with: 

(i) Relaxation effects; 

(ii) Nonlocal effects arising from turbulent diffusion that can give rise to counter-gradient 
transport. 

In virtually all existing full second-order closures for inhomogeneous turbulent flows, 
and pEij are modeled by their homogeneous forms. The pressure-strain correlation is 
modeled as 


$ij - eAij( b) + KMijki{ b) 


dvk 

dxi 


(93) 


as discussed earlier. In Section 5, the equilibrium limit of the Speziale, Sarkar and Gatski 
(SSG) model was provided. For turbulent flows where there are departures from equilibrium, 
the SSG model takes the form (see Speziale, Sarkar and Gatski [2]) 

$ij — —{C\£ + )b{j -f C 2 £ (bikbkj 


where 


--&J kibkiSij) + (C 3 - C;ill /2 )KSij 

— _ 2 — 

-\-C 4 K{bikSjk + bjkSik — -bkiSkibij) 

+CsK{b ik W jk + b jk W ik ) 


(94) 


C x = 3.4, C* = 1.80, C 2 = 4.2, C 3 = - 

0 

c; = 1.30, C 4 = 1.25, C 5 - 0.40, II b = b 

The Launder, Reece and Rodi [20] model is recovered as a special case of the SSG model if 
we set 

Ci = 3.0, C* = 0, C 2 = 0, C 3 = C J = 0, 
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C 4 = 1.75, C 5 = 1.31. 

In most applications, at high. Reynolds numbers, the Kolmogorov assumption of local 
isotropy is typically invoked where 

D^ij ~ ® 

(then, £ij — leSij and a modeled transport equation for the scalar dissipation rate e is solved 
that is of the same general form as that discussed in Section 5). However, this assumption 
is debatable as discussed by Durbin and Speziale [42]. More generally, a representation of 
the form 

D ^ ij ~ 

can be used where the algebraic model (88) of Speziale and Gatski [30], discussed in Section 
5, is implemented. 

The only additional model that is needed for closure in high-Reynolds-number inhomo- 
geneous turbulent flows is a model for the third-order diffusion correlation Cijk- This is 
typically modeled using a gradient transport hypothesis: 

C, jk = V, ]ldmn ^. (95) 

Some examples of commonly used models are as follows: Launder, Reece and Rodi [20] Model 


c^k = -a 


K 


TV 


tm dx m +Tjm 


dr lk drij 

i 1~km ' 


dx. 


dx. 


(96) 


Mellor and Herring [43] Model 

r __ 2 rr R2 ( dT * , dT 'k . dT v) 

' ik 3 ‘ e \ dx, dxj dxk ) 

Daly and Harlow [ 44 ] Model 

C ijk = -2C s -n,pi (98) 

where C s ~ 0.11 is a constant. When these models are used in a full second-order closure, 
counter-gradient transport effects can be described. 

There is no question that, in principle, second-order closures account for more physics. 
This is quite apparent for turbulent flows exhibiting relaxation effects. The return to isotropy 
problem is a prime example where suddenly, at time t = 0, the mean strains in a homogeneous 
turbulence are shut off; the flow then gradually returns to isotropy (i.e., — ■> 0 as t — » oo). 

In Figure 7, results for the Reynolds stress anisotropy tensor obtained from the Speziale, 
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Sarkar and Gatski (SSG) and Launder, Reece and Rodi (LRR) models are compared with 
the experimental data of Choi and Lumley [45] for the return-to-isotropy from plane strain 
(here, r = £ot/Ko). It is clear from these results that the models predict a gradual return 
to isotropy in line with the experimental data. In contrast to these results, all two-equation 
models - including the more sophisticated one based on an anisotropic eddy viscosity derived 
herein - erroneously predict that at r — 0, 6,y abruptly goes to zero. In addition, it is worth 
noting that while the SSG model was derived and calibrated based on near equilibrium two- 
dimensional mean turbulent flows, it performs remarkably well on certain three-dimensional, 
homogeneously strained turbulent flows. The predictions of the SSG and LRR models for 
the normal Reynolds stress anisotropies, compared in Figure 8 with the direct simulations of 
Lee and Reynolds [46] for the axisymmetric expansion, demonstrate this point (here, t * = Tt 
where T is the strain rate). 

While the previous results are encouraging, it must be noted that the Achilles heel of 
second-order closures is wall-bounded turbulent flows: 

(i) Ad hoc wall reflection terms are needed in most pressure-strain models (that depend 
inversely on the distance y from the wall) in order to mask deficient predictions for the 
logarithmic region of a turbulent boundary layer; 

(ii) Near- wall models must typically be introduced that depend on the unit normal to the 
wall — a feature that makes it virtually impossible to systematically integrate second- 
order closures in complex geometries (see So et al. [47]). 

In regard to the first point, it is rather shocking as to what the level of error is in 
many existing second-order closures for the logarithmic region of an equilibrium turbulent 
boundary layer, when no ad hoc wall reflection terms are used. This can be seen in Table 
2 where the predictions of the Launder, Reece and Rodi (LRR), Shih and Lumley (SL), 
Fu, Launder and Tselepidakis (FLT) and SSG models are compared with experimental data 
(Laufer [48]) for the log-layer of turbulent channel flow. Most of the models yield errors 
ranging from 30% to 100%. These models are then typically forced into agreement with the 
experimental data by the addition of ad hoc wall reflection terms that depend inversely on 
the distance from the wall - an alteration that compromises the ability to apply a model in 
complex geometries where the wall distance is not always uniquely defined. Only the SSG 
model yields acceptable results for the log-layer without a wall reflection term. This results 
from two factors: (a) a careful and accurate calibration of homogeneous shear flow (see Table 
3) and (b) the use of a Rotta coefficient |(7i that is not too far removed from one (see Abid 
and Speziale [49]). The significance of these results is demonstrated in Figure 9 where full 
Reynolds stress computations of turbulent channel flow are compared with the experimental 
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data of Laufer [48]. It is clear that the same trends are exhibited in these results as with 
those shown in Table 2 which were obtained by a simplified log-layer analysis. 

The near- wall problem largely arises from the use of homogeneous pressure-strain models 
of the form (93) that are only theoretically justified for near-equilibrium homogeneous tur- 
bulence. Recently, Durbin [35] developed an elliptic relaxation model that accounts for wall 
blocking - and introduces nonlocal effects in the vicinity of walls — eliminating the need for 
ad hoc wall damping functions. While this is a promising new approach, it does not alleviate 
the problems that the commonly used pressure-strain models have in non-equilibrium homo- 
geneous turbulence (the Durbin [35] model collapses to the standard hierarchy of pressure- 
strain models given above in the limit of homogeneous turbulence). The failure of these 
models in non-eq uili bri um homogeneous turbulence can be illustrated by the example shown 
in Figure 10. This constitutes a rapidly distorted homogeneous shear flow that, initially, is 
far fro m equilibrium since SKq/eq = 50 (the equilibrium value of SK/e is approximately 5). 
It is apparent from these 

CHANNEL FLOW 


Equilibrium 

Values 

LRR 

Model 

SL 

Model 

FLT 

Model 

SSG 

Model 

Experimental 

Data 


0.129 

0.079 

0.141 

0.201 

0.22 

&12 

-0.178 

-0.116 

-0.162 

-0.160 

-0.16 

i>22 

-0.101 

-0.082 

-0.099 

-0.127 

-0.15 

t>33 

-0.028 

0.003 

-0.042 

-0.074 

-0.07 

SK/e 

2.80 

4.30 

3.09 

3.12 

3.1 


Table 2. Comparison of the model predictions for the equilibrium values in the log-layer 
(V j e — 1) with the experimental data of Laufer [48] for channel flow. 


HOMOGENEOUS SHEAR FLOW 


Equilibrium 

Values 

LRR 

Model 

SL 

Model 

FLT 

Model 

SSG 

Model 

Experimental 

Data 

bn 

0.152 

0.120 

0.196 

0.218 

0.21 

bi2 

-0.186 

-0.121 

-0.151 

-0.164 

-0.16 

&22 

-0.119 

-0.122 

-0.136 

-0.145 

-0.14 

^33 

-0.033 

0.002 

-0.060 

-0.073 

-0.07 

SK/e 

4.83 

7.44 

5.95 

5.50 

5.0 


Table 3. Comparison of the model predictions for the equilibrium values in homogeneous 
shear flow {V je = 1.8) with the experimental data of Tavoularis and Karnik [50]. 
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results that all of the models perform poorly relative to the DNS of Lee et al. [51]. Even 
the SSG model, which does extremely well for homogeneous shear flow that is not far from 
equilibrium, dramatically overpredicts the growth rate of the turbulent kinetic energy for 
this strongly non-equilibrium test case. 

In the opinion of the author, it is a vacuous exercise to develop more complex models of 
the form (93) using non-equilibrium constraints such as Material Frame- Indifference (MFI) 
in the two-dimensional limit (Speziale [52, 53]) or realizability (Schumann [54] and Lumley 
[16]). While these constraints are a rigorous consequence of the Navier-Stokes equations, they 
typically deal with flow situations that are far from equilibrium (two-dimensional turbulence 
and one or two-component turbulence) where (93) would not be expected to apply in the 
first place. Ristorcelli, Lumley and Abid [55] - following the earlier work by Haworth and 
Pope [56] and Speziale [57, 58] - developed a pressure-strain model of the form (93) that 
satisfies MFI in the 2-D limit. Shih and Lumley [59] attempted to develop models of the 
form (93) that satisfy the strong form of realizability of Schumann [54]. Reynolds [17] has 
attempted to develop models of this form which are consistent with Rapid Distortion Theory 
(RDT). All of these models involve complicated expressions for Mijki that are nonlinear in 
bij. From its definition, Mijki is linear in the energy spectrum tensor Eki(k,t) (see Eq. (47)). 
Since, 

Tfr - iKSn 
2 K 

where T{j is given by (45), it follows that models for Mijki that are nonlinear in bij are also 
nonlinear in Eij . This is a fundamental inconsistency that dooms these models to failure. 
It is clear that is impossible to describe a range of RDT flows - which are linear - with 
nonlinear models (the principle of superposition is violated). Furthermore, Shih and Lumley 
[59, 60] unnecessarily introduce higher degree nonlinearities and non-analyticity to satisfy 
realizability. In the process of doing so, they arrive at a model that is neither realizable nor 
capable of describing even basic turbulent flows (see Speziale, Abid and Durbin [61], Durbin 
and Speziale [62] and Speziale and Gatski [63]). 

Entirely new non- equilibrium models are needed for the pressure-strain correlation and 
the dissipation rate tensor. The former should contain nonlinear strain rate effects and the 
latter should account for the effects of anisotropic dissipation and non-equilibrium vortex 
stretching where Sk 7^ Gk in (84) (see Bernard and Speziale [64] and Speziale and Bernard 
[28]). Models of this type are currently under investigation for the Office of Naval Research 
ARI on Nonequilibrium Turbulence. 
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7. Conclusion 

The following conclusions and recommendations for Naval Hydrodynamics applications 
can now be made: 

(1) For turbulent flows with complex wail bounded or free surface geometries, two-equation 
models with an anisotropic eddy viscosity — that are integrated directly to a solid 
boundary with the no slip condition applied — should be used for the immediate future. 
A new generation of two-equation models, systematically derived from second-order 
closures, has emerged that is far superior to the commonly used K — e model and 
competitive with existing full second-order closures. 

(2) There is no question that full second-order closure models do, in principle, account 
for more turbulence physics than two-equation models. However, current versions of 
these models have major problems when integrated directly to a solid boundary with 
the no-slip condition applied. They also perform poorly in even simple turbulent flows 
that are far from equilibrium. Until these problems are overcome, their use should be 
limited to free turbulent shear flows that are diffusion dominated or to wall bounded 
turbulent shear flows which exhibit complex turbulence physics that does not preclude 
the use of simple law of the wall boundary conditions. 

Research is currently underway, as part of the Office of Naval Research ARI on Nonequi- 
librium Turbulence , to extend these models to turbulent flows that are far from equilibrium 
and to resolve the near- wall problem. With the incorporation of improvements along these 
lines, we should start to see Reynolds stress models make a major impact on the computation 
of the turbulent flows of relevance to Naval Hydrodynamics applications. 
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Figure 2. Time evolution of the turbulent kinetic 
energy in rotating homogeneous shear flow: Com- 
parison of the model predictions with the large-eddy 
simulations of Bardina et al. [37]. (a) 0/5 = 0, (b) 
fi/S = 0.5 and (c) 0/5 = -0.5 (from Gatski and 
Speziale [3]). 
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Figure 3. Comparison of the mean velocity profile in 
rotating channel flow predicted by the new explicit 
ASM of Gatski and Speziale [3] with the experimental 
data of Johnston, Halleen and Lezius [38]. 
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(a) Experiments 



Figure 4. Turbulent secondary flow in a rectangular 
duct: (a) experiments, (b) standard K — e model, and 
(c) nonlinear K — e model of Speziale [22]. 
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Figure 5. Turbulent flow past a backward facing step: 
comparison of the predictions of the nonlinear K — e 
model [22] with experiments, (a) Streamlines and (b) 
turbulent shear stress profiles. 
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Figure 6. Comparison of the predictions of the ex- 
plicit ASM of Gatski and Speziale [3] for skin friction 
with experimental data for the flat plate boundary 
layer (from Speziale and Abid [36]). 
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Figure 7. Time evolution of the anisotropy tensor 
in the return to isotropy problem. Comparison of 
the predictions of the LRR model and SSG model 
with the experiment of Choi and Lumley [45] (from 
Speziale, Sarkar and Gatski [2]). 
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Figure 8. Time evolution of the anisotropy tensor 
in the axisymmetric expansion for eo/TKo = 2.45. 
Comparison of the predictions of the LRR model and 
SSG model with the direct simulations of Lee and 
Reynolds [46]. 
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Figure 9. Comparison of full Reynolds stress calcu- 
lations of channel flow with the experimental data 

(O) of Laufer [48] for Re = 61,600. SSG 

model; FLT model; - • - LRR model; and 

SL model, (a) b n component and (b) 612 com- 
ponent (from Abid and Speziale [49]). 
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Figure 10. Comparison of the SSG, SL and FLT 
model predictions for the time evolution of the turbu- 
lent kinetic energy with the DNS results of Lee, Kim 
and Mom [51] for homogeneous shear flow (SK 0 /e o = 
50). 
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